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ABSTRACT 

We investigate the performance of the parametric Maximum Likehhood compo- 
nent separation method in the context of the CMB B-mode signal detection and 
its characterization by smah-scale CMB suborbital experiments. We consider high- 
resolution (FWHM= 8') balloon-borne and ground-based observatories mapping low 
dust-contrast sky areas of 400 and 1000 square degrees, in three frequency channels, 
150, 250, 410 GHz, and 90, 150, 220 GHz, with sensitivity of order 1 to 10 ^iK per 
beam-size pixel. These are chosen to be representative of some of the proposed, next- 
generation, bolomctric experiments. 

We study the residual foreground contributions left in the recovered CMB maps 
in the pixel and harmonic domain and discuss their impact on a determination of 
the tensor-to-scalar ratio, r. In particular, we find that the residuals derived from 
the simulated data of the considered balloon-borne observatories are sufhciently low 
not to be relevant for the B-mode science. However, the ground-based observatories 
are in need of some external information to permit satisfactory cleaning. We find 
that if such information is indeed available in the latter case, both the ground-based 
and balloon-borne experiments can detect the values of r as low as ^ 0.04 at 95% 
confidence level. The contribution of the foreground residuals to these limits is found 
to be then subdominant and these are driven by the statistical uncertainty due to 
^ , CMB, including E-to-B leakage, and noise. We emphasize that reaching such levels 

>^ ■ will require a sufficient control of the level of systematic effects present in the data. 

1 INTRODUCTION of studies have been performed and published, and which 

have treated the problem on different levels of generality and 

Astrophysical foregrounds are commonly recognized as one detail. The major challenge here is two-fold. Firstly, there 

of the major obstacles on the way to first detecting and later is no general recipe for propagating errors incurred during 

exploiting the scientific potential of the Cosmic Microwave the component separation step, i.e. for including both the 

Background (CMB) polarization signal. This is particularly statistical uncertainty and foreground residual uncertainty 

the case with so called B-mode polarization (Zaldarriaga & Secondly, there is no easily calculable metric measuring the 

Seljak 1997) due to its minute amplitude as compared to impact of the component separation on the B-mode mea- 

the foreground contributions as well as CMB total intensity surement, as both the power spectrum or tensor-to-scalar 

and E-mode polarization signals. In fact current foreground ratio, r, require a proper evaluation of the E-to-B leakage 

models (Page et al. 2007) generally indicate that the fore- (Bunn et al. 2003) . 

ground B-mode signal may be comparable or exceed the ^ucci et al. (2005) have performed a Fisher analysis 
CMB signal by a factor of a few in a broad range of angular ^he problem treating the foreground residuals as an ad- 
scales even in the cleanest available sky areas. Some kind of ditional source of noise, and then estimated the expected 
foreground cleaning or separation procedure will therefore ^ -^^ ^^^^^ approach the starting point was a single 
be necessary and its impact on the final 'cleaned' map of foreground contaminated science channel and a noisy fore- 
the presumed CMB sky needs to be understood and prop- ground template channel from which the level of foreground 
erly taken into account in its subsequent studies. Developing residual was estimated. Ahhough this allows to avoid spec- 
such an understanding is also already of importance for the -^^y^^^ g^.^^^. ^^^^jl ^ foreground cleaning technique, no 
designing and optimization of the future CMB experiments. jjj.gp,. connection exists between the noise values they as- 
This has been recognized for some time and a number sume and properties of any specific experiment. They have 
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also neglected the impact of the E-to-B leakage. A similar 
approach has been followed by Verde et al. (2006), who have 
attempted to link their Fisher matrix considerations to spe- 
cific, fiducial, multi-frequency data sets. The simplified error 
propagation they have adopted implicitly bypasses any real- 
istic component separation approach, and so fails to include 
properly its effects on their final results. They also neglect 
the presence of the E-to-B leakage. Amarie et al. (2005) per- 
formed a Fisher analysis as well, but use specific parameters 
anchored in those of the multi-frequency data set assumed. 
This last work together with Amblard et al. (2007) and Be- 
toule et al. (2009) come the closest in the spirit to what 
we discuss in this paper, although neither of the latter two 
works includes an actual power spectrum estimator account- 
ing for the leakage, what is justified at least in part by their 
focus on full-sky observations. Stivoli et al. (2006) studied 
an application of an Independent Component Analysis based 
approach to the component separation of polarized partial- 
sky maps, resorting to the cleaned map 'pseudo-spectra' as a 
basis for a qualitative assessment of its performance and rel- 
evance for the B-mode work. Dunkley et al. (2009) presented 
a review of most of those earlier approaches, including those 
incorporating a parametric approach similar to the one con- 
sidered in this work, and presented their applications in the 
context of a potential future CMB B-mode satellite mission. 

The approach we propose here is more specific. We fo- 
cus on a particular component separation method and power 
spectrum estimation approach, which wo then use to inves- 
tigate the impact of the foreground separation on the CMB 
B-mode detection and characterization. The component sep- 
aration method is a maximum likelihood (ML) parametric 
approach (Erikscri et al. 2006) in a two-step implementa- 
tion of Stompor et al. (2009). The power spectrum estimator 
is a 'pure' pseudo-spectrum approach introduced by Smith 
(2006) (see also Smith & Zaldarriaga 2007) and elaborated 
on by Grain et al. (2009). Strictly speaking, our results will 
therefore be specific to these two choices. liowcver, given 
that these two methods are working, promising algorithms 
to be implemented in the data analysis pipelines of current 
and future CMB experiments, the results should be of prac- 
tical relevance for many efforts currently going on in the 
field. 

We note also that whenever the firequency scaling laws 
cam be assumed to be nearly perfectly known, as in one 
of the cases we study, and in particular in the small-sky, 
and thus potentially statistics-starved limit, the parametric 
maximum likelihood (ML) method would likely become a 
method of the choice, potentially supplemented by some pri- 
ors, e.g, spatial templates for all or some of the components 
(Efstathiou et al. 2009). The results derived here can there- 
fore be regarded as representative and realistic expectations 
for the performance of classes of the future experiments we 
consider. Moreover, part of the analysis presented here can 
be straightforwardly applied to any component separation 
method in which foreground spectral and amplitude param- 
eters are estimated in separate steps. 

Our focus in this work is on suborbital experiments. 
Those have a potential advantage of selecting the cleanest 
sky areas, but suffer due to the cut-sky effects. They also 
usually have a hmited number of frequency channels with 
which to observe the sky. We consider two kinds of experi- 
ments: those with an access to the high frequencies (^ 250 



GHz) referred to as balloon-borne, and those with access 
limited to frequencies lower than 250 GHz, referred to as 
ground-based. We will also consider some combination and 
extensions of these two cases. We then apply our proposed 
analysis chain to simulated data for different foreground case 
studies, allowing for different levels of mismatch between the 
assumptions made on the analysis and simulation stages, in 
order to evaluate the impact of the component separation 
residuals first on the recovered B-mode power spectrum and 
later on the value of a r which can be derived from such 
data. 

The paper is organized as follows. In Sections 2 and 3, 

we first provide brief descriptions of the specific data anal- 
ysis techniques and their implementation, used throughout 
this paper. In Section 4 we describe our simulated sky model, 
and in Section 5 we define the experimental characteristics 
and a set of foreground case studies. Our results are pre- 
sented in Section 6, and their analysis, concerning the resid- 
uals and their impact on the cosmological B-mode detection, 
is given in Sections 7 and 8, respectively. 



2 PARAMETRIC COMPONENT SEPARATION 
METHOD 

In this section we briefly outline the parametric component 
separation algorithm proposed by Stompor et al. (2009). The 
multi-frequency sky signal is modelled as 

dp — Sp ~\~ Tip, (■^) 

where dp is a vector containing the data from Nfreq fi^e- 
quencies assumed to share a common angular resolution, 
Sp is a vector of Ncomp signal amplitudes to be estimated, 
Ap = Ap {13) is a component 'mixing' or frequency scaling 
matrix with a total of Nspec free 'spectral parameters' /3 also 
to be estimated, and Up is the noise at each pixel p. We can 
write down a likelihood for the data of the form 

-2 In Cdata [s, 13) = CONST + {d-AsfN-^{d-As), (2) 

where N is the noise covariance matrix of the data and we 
have now dropped the pixel index p. This likelihood reaches 
its maximum for the values of s and /3 fulfilling the relations, 

-{A^ffsY N-^ {d- As) =0 (3) 
s = (A* iV-i A) A* iV-i d, (4) 

where ./j denotes a partial derivative with respect to /3i. 
Under the assumption that the spectral parameters are the 
same for a collection of pixels, corresponding to the physical 
assumption that the spectral parameters vary more slowly in 
space than the signal amplitudes, it is possible to substitute 
the generalised least squares solution Equation (4) into the 
likelihood Equation (2), thereby eliminating the sky signals 
s, in order to obtain a spectral index likelihood given by 

-21n£spec (/3) = CONST (5) 
- {A* d)* {A* A)-' {A* d) . 

The spectral parameters that minimize Equation (5) can be 
found using numerical techniques, and then substituted into 
Equation (4) in order to find the corresponding signal am- 
plitudes pixel by pixel. Finally, the noise covariance matrix. 
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describing the properties of the noise contained in the data 
d, is propagated to the component estimates s: 



N, = {A'N-^A)-\ 



(6) 



MlRAMARE^, which is our implementation of this two- 
step algorithm, uses codes from the COSMOMC package 
(Lewis & Bridle 2002) in order to perform an initial conju- 
gate gradient descent to the minimum of the spectral index 
likelihood Equation (5). This is followed by estimation of 
the curvature of the spectral index likelihood on a regular 
grid, which then forms the basis of the 'proposal function' 
for drawing spectral index samples using the Markov Chain 
Monte-Caxlo (MCMC) technique. Once the spectral param- 
eters have been determined from an analysis of the MCMC 
samples, these values are substituted into Equations (4) and 
(6) in order to obtain the signal amplitudes and their covari- 
ance. 

We can check our assumption of the constancy of the 

spectral parameters and the overall 'goodness of fit' by eval- 
uating the log-likelihood Equation (2) at the maximum like- 
lihood value, and comparing it with the number of degrees 
of freedom given by 



r 



(7) 



A poor fit, for instance, due to either the wrong paxametrizar 
tion assumed for the components present in the data or the 

spatially variability of the parameters, will be accompanied 
by an excessive log-likelihood. Other goodness of fit tests are 
also possible and will be discussed elsewhere. 

The ML formalism allows to straightforwardly incor- 
porate the uncertainty due to errors on the calibration of 
the input maps. This can be done by replacing the mixing 
matrix. A, by a product of a diagonal matrix C, represent- 
ing the calibration for each of the single channel maps, and 
the foreground mixing matrix, A. Following Stompor et al. 
(2009), we denote the diagonal elements of C as, oji, so. 



diag(wi)i^o,. 



(8) 



Such a problem is clearly degenerate if no external con- 
straints are imposed on the calibration parameters. In a 
case of actual observations, prior information on the calibra- 
tion uncertainty is in fact usually available and can be often 
expressed as a Gaussian error centered axound some most 
likely value. Mathematically, it just corresponds to multi- 
plying the likelihood in Equation (5) by the relevant priors. 

We note that if the absolute calibration of the result- 
ing component maps is not required one could in principle 
reduce the number of calibration parameters by one by sub- 
suming one of the u>i, say that of the very first channel, 
coo, factors into the overall normalization of the sought-after 
component maps. This is the approach that we will use later 
in this paper, always assuming that the uncertainty of the 
relative calibration of the higher frequency channels with 
respect to the lowest frequency channel can be sufficiently 
well described as a Gaussian random variable with a known 
width. 



3 POLARIZED POWER SPECTRUM 
ESTIMATION 

In this section we briefly describe the approach that we use 
to measure the E- and B-mode power spectra which is based 
on the 'pure' pseudo-C^'s method proposed in Smith (2006). 
Further details of the implementation and tests on simula- 
tions can be found in Smith (2006), Smith & Zaldarriaga 
(2007) and in Grain et al. (2009). The latter work also ex- 
tends the pure approach to the case of cross-spectra. The 
pure pseudo-spectrum estimators retain speed and efficiency 
of the standard pscudo-spcctrum methods and have been de- 
vised to suppress the effect of the E-to-B leakage, thus en- 
suring the near optimality of the estimated B-mode power 
spectrum. The performance of such an estimator is demon- 
strated in the last part of this section where the statistical 
uncertainties on the E- and B-mode reconstruction is eval- 
uated thanks to Monte-Carlo simulation. 



3.1 The pure pseudo-C^'s estimator 

In general a polarization field on the partial sky can be di- 
vided into three classes of modes: pure E-modes, pure B- 
modes and ambiguous modes, which are a mixture of the 
true E and B modes (Bunn et al. 2003). The standard 
pseudo-spectrum approach consists of projecting the polar- 
ization fields 



onto the full-sky harmonic basis of B-modes 



- 



(^-2)! 
{£ + 2)\ 



iid^ - B^) 



(9) 



(10) 



where 8 (8) stands for the spin-raising (lowering) operator. 
However, this basis contains both pure B-modes and am- 
biguous modes on the partial sky. Consequently the stan- 
dard pseudo-C^ estimator includes these ambiguous modes 

in the B-mode power spectra estimates, thus resulting in 
contamination from the much larger E-mode contribution - 
an effect hereafter referred to as E-to-B leakage. This can 
be removed on averaged, but will still lead to a significant 
increase of the estimated spectrum variance. 

The ambiguous modes can however be filtered out by 
projecting the polarization field onto a pure B-mode basis. 
Such a basis is constructed from the spherical harmonics 
and a particular window function W, such that W and dW 
vanish on the boundary of the observed sky: 



-2)! / i(92 _ 
(£ + 2)! 1, d^+B^ 



(11) 



Pseudo-multipoles free of E-to-B leakage can be then com- 
puted from this basis by taking the dot product 



afm = I dn yfj ■ P, 

from which a pseudo-power spectra can be derived 



C 



-B ~B* 



(12) 



(13) 



^ people.sissa.it/~leach/miramare 



The pseudo-power spectrum for the E-modes is identically 
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derived by projecting the Stolces parameters onto the har- 
monic basis for the El-modes, defined as: 

E 1 l ie -2)1 ( d'^ + S' v^^^ 
= 2 V I -i{d' - B') ) (^^^ 

Unbiased estimates for both the E- and B-modc power spec- 
tra are finally obtained by solving the linear system 




where N^, stands for the noise bias and M];.^, for the 
different mode-mode coupling matrices. We emphasize that 
the pure formalism adopted hereafter corrects for the E- 
to-B leakage due to partial sky coverage only, leading to 
elements of M*-"-* much smaller than those of M'^'. How- 
ever, such off-diagonal elements are not strictly zero because 
of pixel-induced E-to-B leakage. This remaining leakage be- 
tween the two type of polarization is nevertheless very small 
and carefully taken into account in the computation of the 
mode-mode coupling matrices. 

The required extensions are included in our implemen- 
tation of the pure pseudo-C^'s estimator, called Xpl're, 
which we use in this work. XPURE is a generalization of the 
XsPECT and Xpol codes (Tristram et al. 2005). It can han- 
dle multiple maps for computing auto- and cross-power spec- 
tra. The different mode-mode couplings due to partial sky 
coverage (i.e. t-to-l' mixing) and pixclization (i.e. residual E- 
to-B leakage) are accounted for during the mode-mode cou- 
pling matrix computation. The code is based on the S^HAT 
library^ - a parallel library allowing for efficient computa- 
tion of spin-weighted spherical harmonic transforms. 

3.2 Sky apodization 

The applicability of the pure pseudo-C^'s estimator strongly 
depends on being able to compute the sky apodizations 
which are needed in the calculation of the relevant pure 

basis functions. Equation (11), and which have to fulfill 
appropriate boundary conditions. Different functions have 
been proposed, ranging from those derived from an opti- 
mization procedure to those based on some analytic expres- 
sions. Their relative merit has been extensively discussed 
by Grain et al. (2009) who showed that the optimized sky 
apodization scheme proposed in Smith & Zaldarriaga (2007) 
leads to the lowest variance on the power spectrum estima- 
tion. 

The underlying strategy for deriving such an optimized 
sky apodization is to search for the W functions which make 
the pure pseudo-C^ as close as possible to optimal, quadratic 
power-spectrum estimator (see Section V of Smith & Zal- 
darriaga 2007). An optimized weighting scheme therefore 
involves a specific sky apodization for each i?-band for which 
the power spectrum is to be estimated, according to the 
signal and noise power aliasing in each band. As shown in 
(Smith & Zaldarriaga 2007), this optimization procedure 
can deal with both external and internal boundaries, due 
to limited sky coverage and holes induced by point source 

2 S^HAT: www.apc.univ-paris7.fr/~radek/s2hat.html 
PUReS^HAT: www.apc.univ-paris7.fr/~radek/pureS2HAT.html 



removal, respectively. Such an optimization procedure how- 
ever assumes that noise and signal arc known. Although this 
is a good assumption for noise and E-modes (which can be 
recovered precisely enough without optimization), this does 
not hold for the B-modes and an erroneous prior may intro- 
duce suboptimality in the B-modes estimate. However, it has 
been shown by Grain et al. (2009) that the B-modes prior 
only mildly affects the performance of the power-spectrum 
estimation using the optimized sky apodization. This weak 
dependence is in fact a direct consequence of the optimiza- 
tion process, which attempts to select an apodization to re- 
duce first the sampling variance due to the El-to-B leakage 
and second the noise variance, with the B-mode variance 
itself usually being subdominant. 

These optimized sky apodizations, leading to the lowest 
statistical uncertainties, will be used throughout this work 
to compute the polarized power spectra from the CMB maps 
estimated from the component separation process. 

3.3 Statistical uncertainties 

A complete characterization of the error budget incurred by 
both the foreground cleaning and power spectrum estima- 
tion processes is mandatory for setting statistically mean- 
ingful constraints on the tensor-to-scalar ratio. We evaluate 
the sampling and noise variance induced by the power spec- 
trum estimation stage using Monte-Carlo (MC) sirrmlation, 
allowing us to determine tlu^ statistical part of the total er- 
ror budget and to demonstrate the performance of the pure 
pseudo-C^'s estimator. 

The observed sky area for the balloon-borne and 
ground-based experiments considered hereafter are shown 
in Figure 2 where holes due to point-source removal are 
carefully taken into account. The sky fraction is roughly 
/sky = 1% and 2.5% for the balloon-borne and ground- 
based experiments, respectively. We assume homogeneous, 
white noise which is deduced from the noise per frequency 
channels using Equation (6). This gives a noise level in the 
CMB maps of approximately 1.65/iK and 3.2/iK per 3.5' 
pixel for the balloon-borne and the ground-based experi- 
ments respectively. The sirrmlations use the WMAP 5-year 
best-fit cosmological model (Dunkley et al. 2009) and the 
simulated B-mode includes the lensing and primordial com- 
ponent with a tensor-to-scalar ratio, r, equal to 0.05. The E- 
and B-mode power spectra are estimated by downweighting 
the simulated maps with the optimized sky apodization as 
described in the previous section and finally binned with a 
band width of A£ = 40, with the lowest £-bin starting at 
I = 20. 

We demonstrate the performance of the pure estimator 
in Figure 1 which shows the input E- and B-mode power 
spectra and the estimated variances of the reconstructed E- 
and B-mode power spectra derived as the standard deviation 
from 500 MC simulations. The three variances displayed on 
this figure arc obtained from, i) the Fisher estimate, or so- 
called /sky-formula, ii) the standard deviation of 500 MC 
simulations using pure pseudo-C^'s estimator and iii) the 
standard deviation of 500 MC simulations using standard 
pseudo-C^'s estimator. 

For the B-modes, the variance is significantly reduced 
by using the pure pseudo-C^'s estimator at angular scales 
where sampling variance is dominating (^ < 700), while the 



Parametric component separation and CMB B-mode detection in suborbital experiments 5 



Balloon-borne experiment 




100 
Multipole t 



1000 



Balloon-borne experiment, r=0.05 



0.1000 r 



a. 



+ 



0.0001 





' — • — ■"'■■'■"■l ' 


J 
















1 ■ 1 ■ ■ - ' ■ 





10 



100 
Multipole t 



1000 



Ground-based experiment 




100 
MulUpole { 



1000 



Ground — based experiment, r=0.05 



0.1000 




0.0001 



100 
MulUpole { 



1000 



Figure 1. Variance of the estimated E-mode (upper panels) and B-mode power spectrum (lower panels) for the balloon-borne (left 
panels) and ground-based (right panels) sky-coverage as shown in Figure 2. From darkest to lightest grey: Fisher estimation (i.e. the 
so-called /gky approximation), standard deviation of 500 Monte-Carlo simulations using the pure pseudo-C^'s estimator as implemented 
in Xpure and, standard deviation of 500 Monte-Carlo simulations using the standard pseudo-C<'s estimator as implemented in XPOL. 
We stress that Xpure provides lower variances than XpOL for the B-modes, while XpOL performs better for the E-modes. 



two approaches lead to similar variances at smaller angu- 
lar scales where noise is the dominant contribution to the 
statistical uncertainties. More specifically, this technique is 
mandatory for the balloon-borne experiment to be able to 
extract the B-mode from the maps of the two Stokes pa- 
rameters at intermediate and large angular scales {£ < 400) 
while the variance is reduced by a factor 2 at those scales for 
the ground-based experiment. Moreover, the pure pseudo- 
Ci's estimator is required for both type of experiments to 
be able to statistically disentangle the inflationary gravita- 
tional waves a.t I < 100 from the secondary, lensing-induced 
B-mode. 



For the E-modes, the standard pseudo-Cf's estimator 
is preferred, for achieving higher accuracy (at least at large 
angular scales for the balloon-borne sky coverage). The rea- 
son for such a higher efflciency of the standard approach as 
compared to the pure one for E-mode is two-fold: on the 
one hand, B-to-E leakage only mildly affects the variance in 
the standard approach and, on the other hand, ambiguous 
modes are mainly composed of E-modes, and so a significant 
amount of information may be lost by removing them using 
the pure approach. 




-44.0} Gtilactio 



Figure 2. Sky areas for the balloon-borne (inner) and ground- 
based (outer, larger patch) experimental setups. The holes mimic 
the effect of masking resolved point sources and are included in 
our optimized apodization calculation. 



4 SIMULATED SKY 

In this section we describe the sky model we use for simu- 
lations in this work. We introduce the 'basic model', which 
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is a simple model of the sky signal consistent with the cur- 
rently available information on diffuse foregrounds, and then 
discuss a number of simple extensions whose physical param- 
eters are poorly constrained by current observations. 



4.1 Sky signals - Basic model 

Our analysis area is centered around RA= 62° and Dcc= 
—45°. This sky region, which is in the anti-sun direc- 
tion during the austral summer, has already been ob- 
served by several CMB polarization experiments including 
Boomerang (Montroy et al. 2006) , ACBAR (Reichardt et al. 
2009), QUAD (Brown et al. 2009) and will be targeted by 
future experiments like EBEX (Oxley et al. 2004; Grainger 
et al. 2008). We first note that WMAP-based estimates of 
the level of unresolved point source power together with 
conservative assumptions about the expected level of radio 
source polarization suggest that this signal will be negligible 
compared to the lensing B-mode signal (Ben Gold, private 
communication). Similar conclusions hold for the infra-red 
sources. Certainly though, a few bright cxtragalactic sources 
will be present in the field both by chance and for the pur- 
poses of in-flight beam mapping and calibration. The process 
of masking out these sources is mimicked in our simulations, 
similar to what done in this respect by Smith & Zaldarriaga 
(2007), as shown in Figure 2. 

Our Galactic sky model on this relatively high Galac- 
tic latitude region consists of two diffuse foreground com- 
ponents: synchrotron and thermal dust emission. The syn- 
chrotron total intensity emission was simulated using the 
408 MHz map of Haslam et al. (1982), with free- free emis- 
sion subtracted and small-scale power added by Giardino 
et al. (2002). We extrapolated this template^ up to 65 GHz, 
using the WMAP five-year maximum entropy method de- 
rived synchrotron spectral index map (Gold et al. 2009) 

The thermal dust total intensity emission is taken from 
the combined COBE-DIRBE and IRAS dust template of 
Schlegel et al. (1998), extrapolated in frequency using the 
scaling accounting for temperature variations as described in 
Finkbeiner et al. (1999). We first extrapolate the dust down 
to 65 GHz for matching our polarization model with WMAP 
data as indicated below, and then extrapolate it back into 
our chosen frequency range according to a uniform greybody 
scaling law inspired by FDS Model 3 (Finkbeiner et al. 1999) 



>ldust OC 



,,/3d + l 



(16) 



where Ta = 18.0K and /3d = 1.65. 

In order to simulate the polarization in this region, we 
normalized the large-scale polarization amplitude to the E 
and B spectra estimated by Page et al. (2007). This nor- 
malization is achieved by first assuming that the polarized 
intensity of both synchrotron and dust is proportional to the 
total intensity, which introduces two free parameters, pdust 
and Payne. lu Order to obtain the large-scale polarization an- 
gles 9, we take the Q and U templates from the WMAP po- 
larized dust template, which is based on information derived 
from starlight polarization and a geometric suppression fac- 
tor taken from a three dimensional Galactic magnetic field 

^ ftp . rssd . esa . int/pub/ synchrotron 



model described in Page et al. (2007). This template intro- 
duces a large-scale modulation to the polarization pattern 
of the synchrotron and dust. We then add Gaussian fluc- 
tuations to the polarization angles on smaller scales follow- 
ing the method of Giardino et al. (2002), assuming a model 

y^cos 20,sin ja ^ o 

Ci ' OC / , where a = — 3. 

With these templates of /, Q and U, we evaluate the 
power spectra for E and B modes in order to finally nor- 
malize the Q and U maps to an effective polarization frac- 
tion. A good match is found for pdust = Psync = 0.1. In 
Figure 3 we show maps of the synchrotron and dust tem- 
plates at 150 GHz. The large-scale modulation introduced 
by the WMAP polarized dust template is clear from the cor- 
related appearance of the synchrotron and dust. A certain 
amount of correlation between these polarized components 
is expected because the Galactic magnetic fleld (GMF) is a 
key ingredient common to both the physics of synchrotron 
emission and dust alignment. Indeed theoretical modelling of 
the GMF is underway by several groups (Page et al. 2007; 
Miville-Deschenes et al. 2008; Waelkens et al. 2009; Jaffe 
et al. 2010). However, a global model of both synchrotron 
and dust polarization in a turbulent GMF at high galac- 
tic latitudes is currently unavailable (see however the recent 
work by Fauvet et al. (2010)), and this is why we adopt 
the semi-empirical approach to our GMF simulations just 
described. Such a modeling can potentially exaggerate the 
overall correlation level between these two components by 
extending it to smaller angular scales. This can in turn have 
important consequences for the performance of the compo- 
nent separation process, as discussed in Section 6.3. 

In the range of frequencies relevant here, the syn- 
chrotron contribution is subdominant compared to the 
dust but increases monotonically with decreasing frequency. 
These two components become comparable at around 70 
GHz, where the total foreground minimum is also found. 



4.2 Sky signals - Extensions 

The problem of insufficient frequency coverage is exacer- 
bated by effects that can generally be described as fore- 
ground complexity, which will increase the biases in the re- 
covered components. We now introduce a few extensions to 
the basic model just described that we will investigate later 
in our analysis. 



4-2.1 Extra small-scale power 

We investigate how increasing the foreground power on small 
angular scales impacts on foreground cleaning by varying 
the parameter a. While in the basic model, the small-scale 

power in polarization is rapidly decaying with a = —3, we 
also study a more extreme case with a = —2. 



4-2.2 Spatially-varying frequency scaling 

Foreground spectral index spatial variations, if poorly esti- 
mated, will compromise the estimation of the CMB signal. In 
turn this compromises the measurement of the CMB power 
spectrum, particularly on the angular scales on which the 
spectral index varies. In the basic model just described the 
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Figure 3. Simulated thermal dust (left) and synchrotron (right panels) Stokes Q and U templates at 150 GHz, for the 'basic model' 
described in Section 4. The high degree of correlation between the two components on large angular scales is imposed by polarization 
angle template we assume. The two foregrounds have roughly the same amplitude around 70 GHz. 
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Figure 4. The dust spectral index variations at 150 GHz adopted 
in this paper (Schlegel et al. 1998). 



dust scales uniformly as Equation (16). We will also investi- 
gate the impact of the FDS temperature variations (Schlegel 
et al. 1998), whose effective powerlaw spectral index (minus 
the mean) at 150 GHz is shown in Figure 4, and whose RMS 
variation is 0.008. 

Approaches for dealing with spectral index variations 
include Taylor expanding the foreground continuum emis- 
sion and fitting the spectral index variations as a new com- 
ponent (Stolyarov et al. 2005), which requires extra fre- 
quency channels. Alternatively, the frequency channels can 
be degraded into many low-resolution, high signal-to-noise 
pixels (Eriksen et al. 2006), in which case an approach for 
dealing with the high-l modes is still required. 

For the multi-pixel approach it is clear that the map 
could be divided up into many sub-regions in which the 
spectral index is estimated, thereby introducing extra fore- 
ground parameters at a cost of increased noise at low £. 
Gaining a quantitative understanding of the impact of spec- 
tral index variations on the multi-pixel method is clearly 
required, since the constancy of the spectral parameters is a 
central assumption in deriving the spectral index likelihood 
Equation (5), but is not the main focus of this work. 



5 MOCK OBSERVATIONS AND 
FOREGROUND CASE STUDIES 

We define here two fiducial experimental setups which we 
will use in the following analysis. They are chosen to re- 
fiect some general characteristics of bolometric experiments, 
but idealized and simplified for demonstration purposes. We 
emphasize that we do not attempt here to forecast a perfor- 
mance of any specific experiment, but rather, on one hand, 
to demonstrate the performance of the considered methods 
in the context of the future B-mode experiments, and on 
the other, to provide reference numbers quantifying the pre- 
cision levels potentially achievable by the small-scale exper- 
iments. 



5.1 Balloon-borne experiment 

We consider a balloon-borne experiment data set with the 
following characteristics: 

• Survey area: Approximately 1% of sky, corresponding 
to around 126, 000 3.5' pixels, showed in Figure 2. 

• Frequency coverage: Three channels at 150, 250, and 
410 GHz each with the same Gaussian beam of FWHM= 8'. 

• Noise level: Homogeneous uncorrelated noise with an 
RMS of 1.5, 4 and 40 fiK on a 3.5' pixel. 

Such a setup can be considered as a minimal choice for a 
balloon-borne observatory. It takes an advantage of not be- 
ing limited in the frequency range by the atmosphere on one 
hand, and on the other it restricts the number of dominant 
foreground components to one. It can be considered as an 
idealization of the data set anticipated from, for instance, 
the EBEX experiment (Oxley et al. 2004; Grainger et al. 
2008). 

We combine these experimental specifications with the 
following three case studies: 

Basic: Based on the basic sky model and the observation 
characteristics as defined above. 

Small-scale power: The basic sky model is augmented 
with extra small-scale power in the dust. 

Varying spectral index: The basic sky model is aug- 
mented with a variation of the spectral index. 
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5.2 Ground-based experiment 

Our ground-based experiment is cliaractcrizod as follows: 

• Survey area: Approximately 2.5% of the sky corre- 
sponding roughly to 320, 000 3.5' pixels, showed in Figure 2. 

• Frequency coverage: Three channels at 90, 150, and 
220 GHz each with FWHM= 8'. 

• Noise level: Homogeneous uncorrelated noise with an 
RMS of 3, 3, and 9 on a 3.5' pixel. 

The frequency range is chosen to fit in the window allowed by 
the atmosphere on one hand, and on the other to conform 
with the optimal working conditions for bolometric detec- 
tors. As a result, the assumed, covered range is quite limited. 
We note that going beyond the lower frequency bound can 
easily be imagined by combining radiometric and bolomet- 
ric data. We will however use the 'minimal' setup defined 
above as our standard case, and extend it on an as-needed 
basis. We finally note that such characteristics are not far 
from those planned for, for instance, the first deployment of 
the POLARBeaR experiment (Lee et al. 2008). 

The qualitative difference between the ground-beised 
and balloon-borne configurations is that the ground-based 
frequency channels are much closer to our foreground min- 
imum at approximately 70 GHz where, by definition, two 
foregrounds are present. This leads us to consider the 
ground-based experimental configurations with a more in- 
volved set of foreground case studies: 

Basic: Based on the basic sky model and the observation 
characteristics as defined above. 

Dust spectral index prior from balloon-borne ex- 
periment: Here we assume the value of the dust spectral 
index determined by the balloon-borne experiment. We dis- 
cuss two cases with all three or only two (150 and 220 GHz) 
channels included. 

Synchrotron template: In this case each channel has 
been corrected for the presence of the synchrotron signal 
via subtraction of a external synchrotron template. We will 
assume that the subtraction is performed down to some pre- 
defined precision level. 

Extra low frequency channel: Here we add an extra 
channel centered at 40 GHz to the basic data set with the 
noise as in the 90GHz map case. 

We point out that only one case above (the basic case) is 
both realistic and self-contained. In all the other cases, the 
presence of some extra external knowledge is postulated. 

In addition to the cases listed above we have also per- 
formed the same analysis assuming a lower noise level of 1/jK 
per 3.5' pixel in all three channels, as well as devised spe- 
cific test runs designed to highlight some particular aspects 
of the performance of the component separation method. 
These include: 

No synchrotron: In this case the problematic syn- 
chrotron emission is removed from the sky model in order to 
understand the way in which it biases the dust estimation 
and subtraction. 

Shuffled synchrotron template: In this case the spa- 
tial morphology of the synchrotron template is randomized. 
The idea here is to understand the effect of accidental corre- 
lations between the foreground components in this regime of 



a restricted number and coverage of the available frequency 
channels. 

We show in Figure 5 the B-modo power spectra of the 
input components at 150 GHz. Depending on angular scale, 
the dust amplitude must be suppressed by between a factor 
of 5 and 25 to have successful measurement of the cosmo- 
logical B-modc signal. 

We note again that the noise levels of the derived CMB 
maps in the basic balloon-borne and ground-based cases as- 
suming a perfect knowledge of the foreground frequency scal- 
ing properties are 1.6 and 3.2 per 3.5' pixel, respectively 
(while in the low noise ground-case we get 0.9 A*K). These 
can be obtained from Equation (6) and are used for the noise 
level shown in Figure 5. 

5.3 Foreground spectral modelling 

Implicit in our approach is the Eissumption that the three 
channel configurations of the ground-based and balloon- 
borne just described will provide information on three, but 
no more than three, parameters: the CMB amplitude and 
the dust amplitude at each pixel, and the dust spectral in- 
dex, as constrained by the ensemble of data and by two 
Stokes parameters. In deriving results in the following sec- 
tion, we begin by assuming exactly the same smooth model, 
Equation (16), for fitting the dust as was used in the simula- 
tions. To some degree this choice is made for expediency, in 
order to define a comparison benchmarks for our case stud- 
ies, and in order to disentangle the cfi^ect of different factors 
on the quality of the CMB estimation. We do however at- 
tempt to gauge the strength of this assumption by studying 
two further case studies that are relevant in this context. 
These are map-level calibration uncertainties and 'spectral 
mismatch'. Intuitively, calibration uncertainties will down- 
grade our ability to infer useful information about the dust 
model and spectral index. Spectral mismatch refers to mul- 
tiplicative factors applied to the dust model at each chan- 
nel, breaking the smoothness of the dust emission spectra. 
This situation could physically occur when the dust scaling 
can not be sufficiently well characterized by the three chan- 
nel experiment. For instance one can conceive of a case in 
which two dust components with warmer and colder tem- 
peratures also have different polarization fractions, leading 
to a sudden change in the dust scaling (and/or position an- 
gle) coincident with the frequency coverage of the exper- 
iment. However, this particular case is thought to be not 
particularly likely for nearby Galactic cirrus, as reviewed by 
Dunkley et al. (2009). Spectral mismatch could also result 
from poorly characterized experimental bandpasses (Church 
et al. 2003). 



6 RESULTS 

6.1 Performance evaluation metric 

To assess the performance of our component separation 
method we will look both at the estimated foreground spec- 
tral indices and at the quality of the recovered CMB maps. 
The latter are clearly not expected to be perfect with po- 
tential contamination arising either due to the noise present 
already in the input, single channel maps, or a failure of the 
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Figure 5. B-mode spectra of the input components at the 150 
GHz channel, for the balloon-borne (upper) and ground-based 
(lower panel) experiments. 



algorithm to perform the separation perfectly. This may re- 
sult either in some level of non-CMB signal still present in 
the map, or in the CMB signal being compromised. These 
two effects are usually referred to as residuals. With the 
noise uncertainty being quite straightforward to character- 
ize using Equation (6), it is our aim to evaluate the level 
of the residuals expected in the foreground case studies and 
then to compare it with statistical uncertainties. The latter 
includes just pixel noise in the case of the maps, and both 
the noise and CMB signal variance on the level of the power 
spectra analysis. 

For each of the case studies described in the previous 
section we first estimate the best-fit spectral parameters by 
maximizing the spectral likelihood. Equation (5), including 
on occasions some extra prior information. Then, given the 
estimated values of the spectral parameters we compute the 
map of residuals as, 



So 



(17) 



where so are the input simulated components and the last 
term on the right hand side subtracts away the noise in the 
recovered components s, given the best-fit values for the 
spectral parameter and the known input noise realization 
n. Directly subtracting the input noise improves upon the 



Table 1. Results for the dust spectral index estimation in the 
balloon-borne case. The last column give the RMS of the total 
foreground contamination left in the CMB map. The input dust 
spectral index is f}^ = 1.65. For comparison, the B-mode signal 
RMS on the pixel scale is 0.19/iK. 



Balloon-borne (/3o 


= 1.65) 




Case 13 


A/3 


RMS (AtK) 


Basic 1.655 


0.009 


0.020 


Small-scale power 1.655 


0.009 


0.021 


Varying spectral index 1.657 


0.011 


0.024 



metric used by Leach et al. (2008) in which residual noise 
was suppressed by smoothing. The CMB element of Equa- 
tion (17), A*~'^^, quantifies the residual foreground signal 
contained in the estimated CMB map, as well as part of the 
genuine CMB signal correctly interpreted by the algorithm 
as the CMB contribution. Finally, we compare the latter 
with the anticipated level of the genuine CMB B-mode signal 
as well as level of its statistical uncertainty due to only CMB 
and noise sampling and cosmic variance. For the last step we 
use as a metric the B-mode power spectra calculated with 
help of the pure estimator. As described in Section 3.3, the 
spectrum variance is estimated via 500 Monte-Carlo simula- 
tions, for which we use the best-fit WMAP 5-year cosmology 
with r — 0.05 as the fiducial model. A satisfactory level of 
foreground cleaning is achieved when the foreground resid- 
ual power spectrum is smaller than the statistical uncertain- 
ties, ensuring that the systematic errors due to foreground 
contamination are subdominant to the global error budget. 

In Section 7 we perform an analysis of these residuals, 
and in Section 8 we express the results in terms of an effec- 
tive detectable value of the tensor-to-scalar ratio, r. 



6.2 Balloon-borne cases 

In Table 1 we report the recovered values of the dust spectral 
index /3 for the three tests we have made: in all the cases 
this parameter was successfully estimated. 

In Figure 6 we show the residuals. A, for the basic case. 
To appreciate the level reached by the cleaning, they can be 
compared with the input foreground contamination at 150 
GHz, shown in Figure 3. The resulting residuals for all the 
balloon-borne cases differ only in a minor way and they are 
always dominated by the unmodelled synchrotron because 
the accurate estimates of the dust spectral index, as derived 
earlier, allow the dust to be subtracted with a superior pre- 
cision. 



6.3 Ground-based cases 

To perform the foreground cleaning, we again assume the 
presence of a single dust foreground because the limited 
number of channels prevents us from performing spectral 
modelling of the synchrotron present in the data. 

The basic result, reported in Table 2, is that the esti- 
mated dust spectral index /3 = 1.875 ± 0.028 is significantly 
biased away from the input value of 1.65, giving rise to resid- 
uals significantly higher than found in the balloon-borne 
cases, as shown in Figure 8 and quantified by an RMS larger 
by a factor ~ 8. This bias of the dust spectral index can be 
explained qualitatively by the fact that in this case not only 
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Table 2. Results for the dust spectral index estimation in the 
ground-based case, as in Table 1. The asterisk indicates the cases 
in which the dust spectral index is fixed, based on the value re- 
covered from the balloon-borne experiment. 



Ground-based (/Jq = 1-65) 



Case 


/9 


A/9 


RMS (/iK) 


Basic 


1.875 


0.028 


0.170 


No synchrotron 


1.642 


0.030 


0.004 


Shuffled synchrotron 


1.667 


0.028 


0.074 


904-150+220 + balloon expt. 


1.655* 




0.077 


150+220 + balloon expt. 


1.655* 




0.034 


External template 


1.682 


0.026 


0.028 



does the synchrotron component, which remains unmodelled 
and unsubtracted, have a higher amplitude than before due 
to a presence of the 90 GHz channel, but also because it 
has a significant spatial correlation with the dust compo- 
nent, as shown in Figure 3 and discussed in Section 4.1. 
We have verified this explanation by making tests first with 
no synchrotron present and later including only the 'shuf- 
fled synchrotron'. The latter case erases the synchrotron- 
dust correlation, artificially converting the synchrotron to a 
white noise-like component with less fluctuations on large 
angular scales. In both these artiflcial test cases, satisfac- 
tory spectral index estimates and foreground cleaning were 
obtained (see second and third panel of Figure 9). More- 
over we found that the pixel by pixel correlation between 
the dust and synchrotron, computed with the Pearson coef- 
flcient C = cov {X, Y) joxOy, has to drop below ~ 15% to 
allow for satisfactory foreground cleaning. 

6.4 Ground-based cases with external information 

The ground-based setup discussed here is therefore not self- 
sufficient and thus unable to provide a appropriately cleaned 
CMB map. In this Section we therefore investigate the ef- 
fect of using 'external information', specifically priors on the 
foreground spectral indices or an external synchrotron tem- 
plate, on the analysis of this data set. 

The first test we made, mimicking a possible real life 
situation, was to impose strictly a dust spectral index prior 
with the value found in the balloon-borne basic case (/3 = 
1.655). At this point, there are no free spectral parameters to 
estimate, and the corresponding least squares components 
can be directly estimated. The residuals for this case are 
shown in the 'delta prior' panels of Figure 10, which demon- 
strates that the dust spectral index prior does help reduce to 
some extent the residuals and the final spectral contamina- 
tion of the B-mode spectrum. Thanks to the high precision 
of the estimation of the dust spectral index in the balloon- 
borne experiment, we find that those residuals are again due 
to the unmodelled synchrotron (on which we will elaborate 
in Section 7). 

Knowing that the 90 GHz channel is contaminated by 
synchrotron, we also have investigated the impact of simply 
dropping this channel, fixing again the dust spectral index 
to the value determined by the balloon-borne experiment. 
Though clearly rather drastic, this choice could in principle 
provide a better foreground cleaning than the three channel 
setup. Unfortunately for the specific case analyzed in this 



work, we found that the remaining two channels arc too 
noisy to produce a CMB cleaned map suitable for any B- 
mode work. We also however have found that two-channel 
setup could be a viable option if the noise in these two chan- 
nels is suppressed down to a ~ 1/iK level - a rather challeng- 
ing goal. Nevertheless this result suggests that some specific 
attention may need be paid to finding the best trade-off be- 
tween frequency bands choice and observation time for this 
kind of experiment. 

We also attacked the problem from the other side of 
the foreground minimum, using information coming from 
lower frequencies. First, we made use of an external tem- 
plate for the synchrotron, whose amplitude is assumed to be 
known with a 10% uncertainty, subtracting it directly from 
the data channels. The satisfactorily foreground cleaned re- 
sult for this case, in which no prior on the dust spectral index 
was assumed, is shown in the second panel of Figure 10. We 
note that though in this test we have assumed a high resolu- 
tion template as only the low-£ modes need to be corrected, 
a low resolution synchrotron templates, as for example an- 
ticipated from the CBASS experiment'' should be sufficient. 
Also, in cases where the overall calibration of the available 
template is considered less reliable than its morphology the 
template marginalization could be a more robust technique 
to be used in this context ( Jaffe et al. 2004) . 

Yet another option we have considered is to extend the 
covered frequency range by adding an extra frequency chan- 
nel operating at 40 GHz. This could be achieved for example 
by co-analyzing the data of two bolometric and radiometric 
experiments observing an overlapping sky area. In our anal- 
ysis the extra 40 GHz was assumed to have a similar noise 
level as the one at 90 GHz. We have indeed found that such a 
combined analysis fares well in terms of residuals amplitude, 
which are comparable to those shown in the second panel of 
Figure 10, in spite of the fact the found best-fit value of the 
spectral index, /? ~ 1.75, is still found to be significantly 
away of the true input value. This fact is illustrated in Fig- 
ure 11 and discussed in Section 7. We note also that the gain 
from the extra channel in the total noise budget of the final 
CMB map. Equation (6), turns out to be negligible. 

To recap the results of this section so far, the balloon- 
borne experiment represents an example of a self-contained 
case where, owing to the fact that its channels are slightly 
displaced from the foreground minimum, a single simplistic 
foreground can be adequately characterized and subtracted. 
It also helps the ground-based experiment to alleviate the 
effect of biases caused by spatial correlation between the 
foregrounds. However it is easy to envisage the need of the 
ground-based experiment for further external information 
on the synchrotron in the form of a template or a lower 
frequency channel. 



6.5 Miscalibration and spectral mismatch 

Here we extend the study conducted so far by incorporat- 
ing two specific systematic effects which are miscalibration 
of the data channels and spectral mismatch between the as- 
sumed and the real model of foregrounds. These two effects 
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Figure 6. Balloon-borne experiment. Maps of residuals in the CMB map, as defined in Equation (17), for the Stokes Q (left) and 
U (right) parameters in the basic case. For comparison, outside the survey border is shown a pure B-mode realization. 
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Figure 7. Balloon-borne experiment. B-mode power spectra of the residuals in the CMB map (cyan curves). As a reference, the 
input and estimated CMB B-mode power spectra are represented by solid black and the solid red curves. The statistical uncertainty of 
the CMB power spectrum estimation is shown by the grey band. For both E- and B-modes and for all these cases, the level of foreground 
residuals is smaller than the statistical uncertainties ensuring a precise enough foreground cleaning. From left to right: i) The basic case, 
ii) Extra small-scale power for the dust emission, iii) Spatially-varying spectral index for the dust emission. 



Table 3. This table summarises how miscalibration errors of the 
two higher frequency channels relative to the lowest frequency 
channel increases the foreground residuals A with respect to the 
basic case. 



Channel miscalibration 


^calibration 




^basic 


0%,2%,2% 


~ 2 


0%,5%,5% 


~ 3 


0%,10%,10% 


^ 4 



differ not only as to their origin, one being due the instru- 
ment properties and the other reflecting our ignorance of the 
physical phenomena relevant to the following case studies. 
They also appear differently within the discussed component 
separation formalism, within which a consistent description 
of only the miscalibration can be incorporated and thus its 
effects properly accounted for. 



Miscalibration 

We simulate a relative calibration error, uncorrelated be- 
tween the channel input maps and applied directly at the 
map level (so that no leakage between the Stokes parame- 
ters occurs) . Though this is clearly a simplification, we note 
that it is not completely unrealistic and may be expected 
for experiments implementing a fast polarization modula- 
tor, such as a continuously rotating half-wave plate used by 
MAXIPOL (Johnson et al. 2007) and under development for 
use with EBEX (Grainger et al. 2008). In such experiments 
the three Stokes parameters can be disentangled from data 
of any single detector, and the resulting maps co-added a 
posteriori in a noise- weighted fashion. The impact of calibra- 
tion errors and uncertainties can be mitigated by modelling 
the calibration parameters at the same time as estimating 
the spectral parameters, as mentioned in Section 2. 

We simulate several cases in which we impose different 
pre-defined miscalibration values, centered on oJi = 1 -I- SuJi, 
introducing Gaussian priors on the calibration parameters, 
centered on uji — 1 with width Suii. We then quantify the 
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Figure 9. Ground-based experiment. B-mode power spectra of the residuals in the CMB map, as in Figure 7. From left to right: i) 
The basic case, ii) A test case with no synchrotron in the simulation, iii) A test case in which the pixels in the synchrotron map were 
reshuffled to remove the spatial correlation with the dust. 



impact of calibration errors by calculating the ratio between 
the residuals in the basic, perfectly calibrated case, and the 
miscalibrated cases. We report these ratios in Table 3 which 
shows the effect of miscalibrating the 250 and the 410 GHz 
channels with respect to the 150 GHz channel, finding for 
instance that 5% calibration errors in these two channels 
leads to foreground residuals that are amplified by a factor 
3. From Figure 7, we can see that this enhancement of the 
foreground spectrum by a factor approximately 10 would 
impact adversely on the large angular scale B-mode estima- 
tion. 



Spectral mismatch 

Here we consider situations where a mismatch between the 
true and postulated scaling laws for the dust component is 
present. In the studied cases, we use different dust scaling 
laws in the simulations, but during the separation process 
we always assume the same, simple dependence as defined 
in Equation (16). The specific laws used in the simulations 
are: the Model 8 from Finkbeiner et al. (1999) and an ar- 
bitrary mismatch, where the dust amplitudes are changed 



Table 4. This table reports how a mismatch applied to the input 
dust models, expressed in percentage per channel, increases the 
foreground residuals with respect to the basic case. 



Spectral mismatch 


^mismatch 
^basic 


0%,0%,3% 


~ 1.5 


0%,3%,0% 


~ 2.5 


0%,0%,7% 


^ 4 


0%,5%,0% 


~ 6 



by some factor from their values as expected in the model 
in Equation (16). Model 8 is a two-temperature model with 
two specific spectral indices and two temperatures for the 
dust. Even if its functional form is different from Model 3 of 
Equation (16), the two models are actually very close in the 
frequency range from 150 to 410 GHz. Fitting the greybody 
scaling to these Model 8 simulations, we found that the final 
foreground residual was small and comparable to the other 
successful cases already shown. 

This motivated us to investigate cases with a larger 
spectral mismatch in which we inserted some discrete mul- 
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Figure 10. Ground-based experiment. B-mode power spectra of the residuals in tiie CMB map, as in Figure 7. Loft: The basic case, 
imposing the the dust spectral index value recovered by the balloon-borne experiment (to be compared with the first panel of Figure 9). 
Right: Subtracting a synchrotron template with an amplitude known to within a 10%. 



tiplicative factors in the dust scaling law used to simulate 
the frequency channels. We progressively broke the assump- 
tion perfect knowledge of the dust spectral behavior, until 
the model is too far from the simulation, leading to B-mode 
biases. In Table 4 we report how the results deteriorate, in 
terms of larger residuals A*"^^^ , for some mismatch choices 
in different channels. The basic result is that mismatches of 
upwards of 5% in the dust scaling can lead to an enhance- 
ment by a factor 6 and upwards of the foreground residual 
level. It is perhaps not surprising that effective modelling 
and subtraction of foregrounds using few channels and few 
free parameters depends on the underlying smoothness of 
the foreground frequency scaling. 

We note that this test differs from the miscalibration 
case discussed earlier, as only one of the components am- 
plitude is modified and no prior is used in the separation 
process. 



7 ANALYSIS OF THE RESIDUALS 

In the simulation environment we have the power to con- 
trol details of all the aspects of the separation process. In 
this Section we take the advantage of this fact and inves- 
tigate the nature and origin of the residuals A shown so 
far in the paper. We emphasize that the considerations pre- 
sented below do not depend on how the estimate of the 
spectral parameters have been obtained, and therefore they 
apply more generally than just to the specific parametric 
ML estimator considered here. In fact, the analogous rea- 
soning could be applied, and similar conclusions drawn, in 
a case of any two step approach in which first the spec- 
tral indices estimates are derived and then the sky com- 
ponents estimated via Equation (4). This includes FastICA 
(see Bottino et al. 2010, and reference therein), neural net- 
works (N0rgaard-Nielsen & j0rgensen 2008), and Correlated 
Component Analysis (Bonaldi et al. 2007). 

As introduced in Equation (4), the operator we apply 
to the data set d to recover the component estimates s is, 

W {(3) = (A* (/3) N-^ A{I3))-' A' {(3) N'' , (18) 

where we explicitly highlight the dependence on the recov- 
ered dust spectral index /3. Neglecting the presence of the 
noise, the data can be written as d = A{f3o)so (hereafter the 



subscript refers to true, i.e., input rather than estimated 
quantities), and therefore, 

s = Wd ^WA{(3o) so = Z{l3) so, (19) 

and thus the residuals can be written down as, 

A = s -so = (Z(/3)-/) So, (20) 

where the last definition of A coincides with Equation (17) if 
no noise is considered. Here So refers to the true, input com- 
ponent, which is modelled in the separation process and is 
thus a subvector of so- The matrix / is made of a square unit 
matrix, corresponding to all modelled components, supple- 
mented by extra columns of zeros - one for each unmodelled 
component. 

The size of the matrix, Z, depend on the number of 
actual and derived sky components and not on the number 
of the observed frequency channels. In most of the examples 
shown in this paper, Z' is a 2 x 3 matrix, since we attempt 
to recover only two (out of three) components. We note that 
once we have estimated /3, the operator Z {(3), that trans- 
forms the input components into the output ones, can be 
readily computed since, in the simulations, we know the in- 
put scalings. In this case. Equation (20) provides insight into 
the origin of the residuals and their relative amplitudes. 

We first observe that, 

Z (/3 = A.) = 1, (21) 

if the number of assumed and actual components is the 
same. If there are more components used for the simulations 
then subsequently recovered, this will no longer be the case. 
However even then the maximal square block of the matrix 
Z {(3 = (3o) will be equal to a unit matrix. This is shown in 
the upper part of Table 5. Clearly, even perfect knowledge of 
the true dust spectral index does not assure the lack of the 
residual. Nevertheless, in such a case each of the recovered 
components contains only a contribution of this component 
plus one due to the unmodelled signal. In a specific case 
of the balloon-borne experiment considered here nearly all 
of the unmodelled synchrotron is added mostly to the re- 
covered CMB signal, given the similar scaling of both these 
components in the respective frequency bands. 

The matrix Z{P) computed in a more general and real- 
istic case, when the spectral index is unknown and needs to 
be estimated from the data, is shown in the middle part of 
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Table 5. We note first that as before, and for tlie same rea- 
son, the unmodelled synchrotron contributes predominantly 
to the CMB residuals. Nevertheless, the dust is now divided 
in between the two recovered components, though the dust 
signal found in the recovered CMB is very subdominant. 
Also in the previous case, the recovered CMB component 
contains the entire CMB signal, which is completely absent 
in the recovered dust template. This is reminiscent of Equa- 
tion (21), which in the present case holds only for the CMB 
component reflecting the fact that the perfect black-body 
derivative scaling is assumed on both the simulation and re- 
covery stages. In this case again, the CMB residuals, A'"^^^, 
do not contain any CMB. This is no longer true if the cal- 
ibration errors are allowed for as shown in the bottom of 
Table 5. In this case the recovered CMB component con- 
tains only part of the total CMB signal as determined by 
W-matrix weighted average of the relative calibration errors 
for each of the channels. The remainder of the CMB is then 
found in the recovered dust. We note however that though 
for the calibration errors considered here these effects are 
small, the recovered CMB residual is now indeed typically a 
mixture of the CMB, dust and synchrotron signals. We point 
out that, maybe somewhat counterintuitively, the elements 
of Z in any column do not have to sum up to unity. This 
reflects the fact that due to our wrong assumptions about 
the spectral parameters values, the estimated components 
contain overall 'more' of the input components than there 
really is. This effect is small, if the estimated /3 values are 
close to the true ones, as to first order in /3, the columns of 
Z do sum up to (nearly) unity. 

Finally, in all the cases we find that although the code 
leaves in most of the unmodelled synchrotron, it cleans the 
dust to better than ~ 0.5%. The latter depends on the spe- 
cific value of j3 assumed for the recovery as shown in Fig- 
ure 11, which shows the relative contamination of the recov- 
ered CMB due to the dust as a function of the recovered 
spectral index in the case of the balloon-borne and ground- 
based basic cases. The shaded band indicates the precision 
which is needed to avoid contamination in the cosmological 
B-mode recovery. 

The results obtained in the ground-based cases are qual- 
itatively similar to the ones described above. However, in the 
basic case, due to the different assumed frequency coverage, 
even for the true value of /3 we find non-zero contributions 
of the synchrotron in both recovered components. Moreover, 
the contribution to the CMB is more than three times that 
of the actual synchrotron signal at 150 GHz. Assuming in 
turn the best-fit value (5, we find that the synchrotron levels 
in both recovered components remain essentially unchanged, 
however the absolute dust contribution to the CMB template 
increases to become of the same order as that of synchrotron. 



8 COSMOLOGICAL B-MODE DETECTION 

The central question of this work is to understand how much 
the recovery of the cosmological B-modes is affected by the 
presence of the foregrounds or foreground residuals left over 
by some foreground cleaning technique. This question is of- 
ten phrased as a question about the detectable values of 
the tensor-to-scalar ratio, r, of the corresponding primor- 
dial power spectra, and the calculated limiting values of r 




Figure 11. Dust suppression factor in the 150 GHz channel of 
the basic case of the the balloon-borne (solid line) and ground 
based (dashed line) experiments, as a function of the recovered 
dust spectral index f}^. The horizontal shaded band is indicative 
of the requirement on the dust suppression factor in order to do 
B-mode science. 



Table 5. This table reports how the input components are 
weighted in the outputs calculated using Z as defined in Equa- 
tion (19), for the balloon-borne basic case and an ideal case where 
the dust spectral index is known. The RMS values of the CMB, 
dust and synchrotron signal in the studied patch are 3.1, 0.6, 0.02 
/iK, for CMB, dust and synchrotron respectively. 



Balloon-borne: Ideal case, 


l3 = l3o = 1.65 


Input; 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


1.000 


0.000 


1.003 


Dust 


0.000 


1.000 


-0.037 


Basic case, /3 = 1.655 


Input; 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


1.000 


0.005 


1.003 


Dust 


0.000 


0.994 


-0.036 


5% Miscalibration case 


/3 = 1.639 


Input; 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


0.988 


0.005 


0.992 


Dust 


-0.0004 


0.982 


-0.037 



are dependent on the implicit and explicit assumptions made 
in the course of its derivation. 

In the context of this paper this question is, however, 
well defined as we study specific foreground separation and 
power spectrum estimation algorithms. Our goal here is to 
derive values of r, which can not only be detected by the con- 
sidered experiments but also convincingly argued for from a 
perspective of an observer as indeed being driven by the pri- 
mordial signal of the cosmological origin. The limits we aim 
for here are not to be considered as some 'ultimate' lower 
limits on r, as often quoted in the literature, but rather rep- 
resentative of the potential of the specific experiments and 
data analysis techniques considered here. Reaching values of 
r lower than our estimates by the actual experiments once 
they are deployed and operating, may not only be plausible 
but indeed is expected owing to the build up, over time, of 
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Table 6. This table reports how the input components are 
weighted in the outputs, for the studied ground-based cases. The 
RMS values of the CMB, dust and synchrotron signal in the stud- 
ied patch are 3.1, 1.3, 0.02 /iK, for CMB, dust and synchrotron 
respectively. 



Ground-based: Ideal case, 


/3 = /3o = 1.65 


Input: 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


1.000 


0.000 


3.116 


Dust 


0.000 


1.000 


-1.45 




Basic 


case, P = 1.875 


Input: 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


1.000 


0.067 


2.993 


Dust 


0.000 


0.919 


-1.450 


Synchrotron template, /3 = 1.682 


Input: 


CMB 


Dust 


Synchrotron 


Output: 








CMB 


1.000 


0.004 


3.108 


Dust 


0.000 


0.995 


-1.441 



our knowledge of foreground properties and data analysis 
tools. 

We start by developing a model to describe statistically 
the foreground residuals found in the recovered CMB map, 
thus extending the pixel-domain discussion of the last sec- 
tion into the power spectrum domain. We assume that our 
estimates of /3 are unbiased and that the obtained statistical 
uncertainty of the (3 recovery is small, as is indeed the case 
in the cases considered. Now, if all the actual components 
are included in our data model, as defined by W (/3), it is 
accurate to write, 

,dZ{l5o) 



Z(/3). 



5/3- 



(22) 



where 5/3(s /3 — /3o) is assumed to be a Gaussian variable 
with a zero mean and the dispersion as derived earlier, and 
we used the fact that Z (/3o) = 1, as in Equation (21). Using 
Equation (20) we can then express the residuals of all the 
modelled components as, 



. mod 



[Z{I3)-1] so 
dW{/3o) 



5f3 



dp 



A{l3o) So, 



(23) 



(24) 



and use a first row of this matrix equation to compute the 
foreground residuals as found in the recovered CMB map. 
We first introduce a tensor, , defined as, 



df3k 



(25) 



and we can then express a combined residual in the CMB 
map {i = 0) due to all the modelled, non-CMB component 
as, 

A^^^ = Y,5f3,c.l^4. (26) 

k,j 

This shows that the residuals behave like templates with 
amplitudes randomized due to the impact of the CMB and 
noise variance on the determination of the spectral param- 
eters. Hereafter we will assume that the latter are Gaussian 
variables centered at the true value, /3o, and with dispersions 



as estimated from the data, S'^(= ((5/3(5/3')). Consequently, 
for pure spectra averaged over the statistical ensemble (of 
the noise and CMB) we can then write. 



k.k' j.j' 



(27) 



where Cq"'^ are the (pure) auto- and cross- spectra for every 
pair of the actual j and j' non-CMB components. 

In a case of foreground components which can not be, 
and/or are not modelled in the separation process, we can 
no longer use the procedure outlined above to estimate their 
residuals. From the previous section we note that such resid- 
uals depend only weakly on the assumed spectral parameters 
and thus will not be stochastic over the CMB plus noise re- 
alizations and should be rather treated as a bias. Moreover, 
in the specific cases considered in this paper we note that 
the element of the matrix, Z, which determines the weight 
with which the component is added to the cleaned CMB 
map is on order of at most a few, what together with the 
fact that the overall expected level of the synchrotron is very 
low, leads to the conclusion that indeed it can be neglected 
for the estimations of r as derived here. This last statement 
can be phrased somewhat more formally by using the ap- 
proach of Amblard et al. (2007), which allows to quantify 
values of r, denoted as Vres, which will not be affected by a 
bias due to some present residual. Indeed, we define the two 
quantities: 



s(r) 



■cmb 
I 



ir), 



(28) 
(29) 



In the successful cases considered in the paper we find that 
typically u ~ s{rres) for r^el = 0.005 and ril], = 0.015 
for the balloon-borne and ground-based experiments respec- 
tively. They represent the smallest values of r for which the 
unmodelled residuals due to synchrotron can be neglected. 

We will now assume that the bias due to the unmodelled 
components is negligible and proceed to the estimation of 
the parameter r accounting for the extra variance due to 
the modelled component residuals. For this we use a Fisher 
like approach, which in our case reads. 



F{r 



dCi: 



dr 



dr 



(30) 



where b denotes the £-hin number used for the power spec- 
trum estimation, and. 



■'bb' 



(r) = Var(cSr«+"°'='^) +Cfe^^^+°°'='^Cb^ 

, >--tCMB+noise /^A , >--yA 



(31) 



Here the first term is the covariance of estimated pure CMB 
B-mode spectra, which we assume to be diagonal in the 
bin-space and which is estimated via MC simulations for 
a grid of values of r. (^CMB-i-noise estimated pure 

CMB-|-noise spectra also computed on a grid of r values. The 
partial derivatives are computed using the binned theoreti- 
cal spectra, obtained in this case from the camb code. This 
is justified given that the pure estimator is assumed to be 
unbiased. The last term in Equation (31) describes the vari- 
ance due to the residual foreground treated here as a tem- 
plate fully correlated between different bins. The third and 
fourth terms reflect the cross-terms between the foreground 
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residuals and CMB+noise. These again are non-diagonal in 
the bin domain. We note that the off-diagonal template-like 
correction appearing in Equation (31) corresponds roughly 
to the tensor-product term found in Equation (A4) of Stom- 
por et al. (2009), which indeed describes the correction to 
the noise correlation matrix, Equation (6), due to the fore- 
ground residuals. 

We then search for the values of r^j such that r > rd — 
2 F~^^^ (r), which can hence be detected at a confidence 
level not smaller than 95%. For the experimental setups de- 
scribed in this work, we find that ~ 0.04 both for the 
balloon-borne and the ground-based cases. This value of r 
is sufficiently high that the bias expected due to the unmod- 
elled synchrotron is indeed much smaller than the detected 
values. We also find that the derived value of r is limited by 
the CMB-(-noise variance, with the foregrounds effects be- 
ing subdominant. This observation is perhaps not surprising 
given the low level of the residuals as discussed in the previ- 
ous Sections. However, only an analysis as the one presented 
in this Section, can properly account for the bin-bin corre- 
lations of the foreground residual templates. 

We note that up to this point we have assumed that 
we are privy to some insights as to the true nature and 
morphology of the foreground signals and their scaling well 
beyond what is usually available to an actual CMB observer 
and data analyst. In a real life situation we will lack some 
of that information. Specifically we will be likely ignorant of 
the true value of the /3 parameters, i.e., /3o, and cross-spectra 
of all the actual sky components, C^'-'^ . That may not look 
like a big issue given our conclusion above stating that the 
dominant uncertainty will be due to the CMB and noise 
variance. However, this conclusion may need to be corrobo- 
rated case-by-case using data, analyzed together with some 
necessary external information. Whenever it turns out not 
to hold then a consistent procedure to account for the extra 
effects may be needed. 

In the case at hand this can be done by replacing the 
needed information with their best estimates derived in the 
analysis process. We will thus use the best-fit value in place 
of (3o and the pure spectra of the components estimated 
as a result of the separation process, and corrected for the 
noise bias using Equation (6) in the case of the auto-spectra, 
j = j' . In Figure 12 we show how this procedure fares in 
the case of the dust component in one of the ground-based 
case example considered in this paper. To evaluate the bias 
due to the unmodelled synchrotron, we need to use external 
data. We could correct for such a bias, if the available data 
are of sufficiently high precision, and include the resulting 
variance in the final power spectrum error budget. However, 
more typically, one would rather aim to show that the bias is 
indeed negligible in the sense defined above in Equation (29). 
For this task the required external data however need not 
being very precise, not least due the fact that the unmod- 
elled components considered here are assumed to be truly 
subdominant. 

Using the estimated quantities in our r-estimation pro- 
cedure we recover the limits on r essentially identical to 
those found before. We thus conclude that values of r ^ 0.04 
are not only detectable at greater than the 95% confidence 
level, but can be detected in the realm of an actual experi- 
ment and argued for as of a cosmological origin, all of that 
providing a sufficient control of systematic effects. 



Residuals 




20 100 1000 

Multipole 



Figure 12. Comparison between actual and estimated residu- 
als for the ground-based, no synchrotron case, calculated using 
Equation (27). 

9 CONCLUSIONS 

In this paper we study the performance of the maximum 
likelihood parametric component separation method from 
the point of view of its application to the CMB B-mode 
polarization analysis. We investigate the residuals left over 
from the separation in both the pixel and harmonic domains. 
We propose an efficient framework for evaluating the pixel 
domain residuals in the simulation, and show how it can 
be used to gain important insights into the separation pro- 
cess. We then compute the power spectra of the recovered 
CMB maps, as well as maps of the residuals, using the pure 
pseudo spectra technique, and estimate their variances using 
Monte Carlo simulations. Finally, we propose a Fisher-like 
approach to evaluate the effects of the foreground residuals 
on the r parameter and use the latter to derive some esti- 
mates of typical values of r, which are potentially detectable 
by the considered experiments at the 95% confidence level. 
The latter estimates thus include the uncertainties due to 
sampling variance, noise scatter, E-to-B leakage, and fore- 
ground residuals, all of which are consistently propagated 
through the proposed pipeline. 

We focus here on the small-scale, bolometric experi- 
ments, broadly dividing them into two classes, referred to 
as balloon-borne and ground-based setups, both observing 
in three different frequency bands. We find that the balloon- 
borne case, with frequency bands at 150, 250, and 410 GHz, 
provides a robust experimental setup for the detection of 
the B-mode polarization. The foreground residuals in the 
recovered CMB maps derived in this case are found to be 
usually subdominant. This is true whenever the assumed 
data model is indeed correct, but it also holds when some 
small systematic effects are permitted. Selected effects of 
this kind considered in this work include relative calibration 
errors, unmodelled spatial variation of the spectral parame- 
ters, and spectral mismatch between assumed and true spec- 
tral scaling laws. We emphasize that all these systematics, 
though manageable if sufficiently small, may lead to spuri- 
ous effects in general, and therefore need to be controlled in 
actual experiments with a sufficient precision, which need to 
be determined specifically for any experiment. 

The success of the considered balloon-borne cases is 
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related to the wide frequency range available to such ex- 
periments, which permits selecting frequency bands at the 
sufficiently high frequencies to avoid the unwanted residual 
synchrotron. The latter is found to be a dominant source of 
the bias for the ground-based experiments. In the balloon 
case we can also afford a long leverage arm between the 
lowest (CMB-dominated) and the highest (dust-dominated) 
frequency bands, which plays a pivotal role in setting tight 
constraints on the spectral parameters of the dust. From our 
Fishcr-likc analysis, wc show that one could detect r values 
as low as 0.04 at the 95% confidence level with such experi- 
ments, if both our models and measurements are sufficiently 
well characterized. 

For the ground-based case the atmospheric loading re- 
stricting the available frequency window proves to be a sig- 
nificant limitation. We find that even in an absence of any 
systematic effects with three frequency bands set at 90, 150, 
and 220 GHz, it is generally not possible to produce suf- 
ficiently clean CMB maps. This is due to the unmodelled, 
and thus not separated, synchrotron contribution, which is 
significant enough (if the polarized emission is at the level 
suggested by WMAP) at these frequencies to bias the esti- 
mation of the dust spectral parameters. Wc point out that 
this contribution has been neglected in some earlier works, 
which consequently has arrived at a different conclusion. 
This therefore emphasizes the importance of accurate sky 
modelling for this kind of the analysis. Nevertheless, we find 
a satisfactory cleaning can be achieved in such a case if 
some external information is available. In particular, we dis- 
cuss the extended ground experiment analysis allowing for 
the presence of the extra lower frequency channels, rough 
synchrotron templates, and priors on the dust scalings. We 
find that in such cases, and under realistic assumptions, the 
ground-based experiments should reach a sensitivity roughly 
matching those found in the balloon-borne case, in terms of 
a detectable r parameter. We also conclude that, for both 
these types of the experiments, the foreground contamina- 
tion anticipated in low-contrast dust regions, should not be 
an obstacle preventing them from exploring the parameter 
space of r down to the values ~ 0.04. Indeed, for the consid- 
ered experimental setups this limit is determined by the un- 
certainty due to the CMB itself and the instrumental noise, 
with the effects of the residual foregrounds found to be sub- 
dominant. In the realm of the actual observations, whether 
these limits are reached will be crucially dependent on the 
control of systematic effects. Wc note here however that the 
limits derived on r arc strongly dependent on the level of the 
noise assumed in the input single-channel maps. These can 
be therefore improved upon, if a deeper integration of the 
same field is performed. However, if no additional external 
information is used, those limits will remain appropriately 
higher than the Tres values obtained earlier, and below which 
foreground bias would become significant. 

In this context we point out that the framework de- 
scribed in this paper provides a blueprint for similar studies 
focused this time on systematic effects. It can be also ex- 
tended to perform a realistic experiment optimization pro- 
cedure from the viewpoint of detection of the B-mode signal 
of cosmological origin. 
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